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ABSTRACT 

This paper concerns the image reconstruction from a few 
projections in Computed Tomography (CT). The main objec- 
tive of this paper is to show that the problem is so ill posed 
that no classical method, such as analytical methods based 
on inverse Radon transform, nor the algebraic methods such 
as Least squares (LS) or regularization theory can give sat- 
isfactory result. As an example, we consider in detail the 
case of image reconstruction from two horizontal and verti- 
cal projections. We then show how a particular composite 
Markov modeling and the Bayesian estimation framework 
can possibly propose satisfactory solutions to the problem. 
For demonstration and educational purpose a set of Matlab 
programs are given for a live presentation of the results. 

1. INTRODUCTION 

The image reconstruction problem in Computed Tomogra- 
phy (CT) is presented in Fig. 1 and Fig. 2 shows the forward 
modeling of the problem via the Radon Transform (RT) and 
the basics of the analytical RT inversion based methods such 
as Filtered backprojection ClIlEllllEillTlIHl. 
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g^{rur2) = / f{x,y,z) dl g^{r) = / f{x,y) dl 

Forward problem: f{x,y) ox f{x,y,z) — > g,^{r) or g^{ri,r2) 
Inverse problem: g^(r) org(^(ri,r2) — > f(x,y) or f{x,y,z) 

Fig. 1 : Tomography X 



However, it is so evident that these methods cannot give 
satisfactory results in cases of very limited number of pro- 
jections as is the case we consider in this papaer. To be 
able to introduce the necessary prior information needed to 
overcome the lack of information in the data, we consider 
the algebraic methods. Fig. 3 shows the discretization step 
of the forward problem which transforms the linear continu- 
ous RT equation to a system of finite linear equations which 
is g = Hf. It is then evident that this system is under- 
determined and that the problem has an infinite number of 



solutions. As a demonstrative example, we consider the 
case of image reconstruction from only two projections and 
study the structure of the matrix H in this particular case 
and show easily that neither the neither the minimum norme 
least squares (MNLS) nor the generalized inversion and nor 
the quadratic regularization |9| can give satisfactory result to 
this problem. 
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Fig. 2 : X ray Tomography and Radon Transform 



We also show that even applying the positivity constraint 
is not enough to obtain satisfactory results and that there is 
a need for more informative prior knowledge. Finally, as the 
main contribution of this work, we show that the Bayesian in- 
ference framework and the composite Markov modeling can 
possibly be of great help to develop new reconstruction meth- 
ods whith possibly satisfactoy results. We consider in par- 
ticular a composite and hierarchical Intensity-labels Markov 
modeling with a Gauss-Markov modeling for the intensity 



field and a hidden Pottz Markov field for the region labels and 
propose new reconstruction methods which can be applied in 
many imaging systems, and particularly, in Non Destructive 
Testing (NDT) imaging applications. 

2. DISCRETIZATION OF THE PROBLEM 

As we mentionned before, for demonstration purpose, we 
conside the particular case of image reconstruction from only 
two projections. Also, for the sake of simplicity, we give de- 
tails about a very reduced case of a (4 x 4) pixels image. 
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Fig. 3 : Discretisation du probleme 



If we note by the vector / = [/i , • • • ,/i6]^ the pixels val- 
ues of the image f{x^y) and by the vector g = [gi, • • • jgs]^ 
the values of its two projections g{r^^) along the horizontal 
= and vertical = 90 angels, and assuming Ax = 1 , = 
1 , Ar = 1 , the we have: 
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Al=[ 1111000000000000; 
0000111100000000; 
0000000011110000; 
0000000000001111]; 

A2=[ 0001000100010001; 
0010001000100010; 
0100010001000100; 
1000100010001000]; 

A= [Al; A2] ; 

f = [ ; 1 1 ; 1 1 ; ] ; 

P=A^f (:) ; 
gives : 

g' = [02200220] 

Thus, we have modeled the forward problem. Now, we are 
going to consider the inverse problem which is given g find 

/■ 

It is evident that this inverse problem is under- 
determined, and that it has an infinite number of solutions. 
Here are four examples: 
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3. BACKPROJECTION AS THE ADJOINT 
OPERATOR OF RT AND ITS EQUIVALENT MATRIX 
TRANSPOSITION 

Comparing the continuous RT and its corresponding dis- 
cretization and the adjoint Backprojection operator given in 
Fig. 2, and its corresponding discretization, we see easily 
the equivalence of Backprojection and the transposition of 
the the matrix A. Thus, / = A^g corresponds to an image 
obtained by backprojecting the projections and it is easy to 
see that: A^ = [ A\ \ A2 ] Thus, computing this solution is 
easy: 

fh=A' ★p; reshape (fh, 4, 4) 



Noting also by gi = [gi,-",^4]^ = [giw ' ' ^guY , 9i = 
[gsr-- ,gsy = [giw" ^gi^y and the matrices Ai, A2 and 
A such that: 



gi=Aif, g2 = A2f, g = Af 



Then, considering the image 
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We may also note that / 
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the following Matlab code lines: 



each being the backprojection of each projection. We may 
also remark that this solution is not very far from the result 



of the convolution of the original image with the following 
impluse response 

1 

1 2 1 
1 

4. LS, MNLS AND GENERALIZED INVERSION 

Let consider the two symetric matrices A and AA^: 
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4000 1 1 1 1 
0400 1111 
0040 1 1 1 1 
0004 1111 
11114000 
11110400 
1 1 1 1 0040 
1 1 1 10004 
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A'A- 



2 1 1 1 1 000 1 000 1 000 
12 110 1000 1000 100 
1 1 2 1 00 1 000 1 000 1 
1 1 1 2000 1 000 1 000 1 
1000211110001000 
0100121 101000100 
0010112 10010001 
000 1 1 1 1 2000 1 000 1 
1000 10002 1111000 
0100010012110100 
0010001011210010 
000 1000 1 1 1 12000 1 
1000 1000 10002 1 1 1 
1 000 1 000 1 00 1 2 1 1 
0010001000 101121 
0001000100011112 



and compute their singular values: 

AAt=A^A' ; svd ( AAt ) 

svd(AAO = [84444440] 

AtA=A' ^A; svd (AtA) ; 

svd(A'A) = [8444444000000000] 

We can then remark that both are singular. 

We may remind that the least squares (LS) solutions are 
defined as 

f = argnnn{\\g-Aff}, 
and if A^A was invertible, then we had: / = {A^ A)~^A^g. 



as: 



In the same way, the minimum norme solution is defined 
/ = arg min {\\ff} 



and if AA^ was invertible, then we had: / = A^ {AA^)~^g. 

As we noticed, the two matrices AA^ and A^ A are sin- 
gular and thus we can not define those solutions. However, if 
we only consider their diagonal elements, we can define: 

fh=diag (1 . /diag (AtA) ) ★A' ★p; reshape (fh, 4, 4) 
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fh=A' ^diag (1 . /diag (AAt) ) ★p; reshape (fh, 4, 4) 
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Also, we can use the technic of truncation of the singular 
values to define an unique generalized inverse solution 



<g,Uk > 

L — ^ — 



where and are, respectively, the eigenvectors of AA^ 
and A^ A and their corresponding eigen values. 
[U, S, V] =svd (A) ; 

s=diag (S) ; sl= [ 1 . /s ( 1 : 7 ) ; zeros (1,1)]; 
Sl= [diag (si) ; zeros (8,8)]; 
fh=V^Sl^U' ★p; reshape (fh, 4, 4) 

In this example, K = 1 and the GI solution can be com- 
puted by 

f h=svdpca (A, p, . 1 , 7 ) ; reshape ( f h, 4 , 4 ) 



/ = 
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0.2500 0.7500 0.7500 
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-0.2500 0.2500 0.2500 
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or by the following iterative algorithm: 

for 1^=1:100; 

fh=fh+. l^A' ^ (p-A^fh ( : ) ) ; 
end; 

reshape ( f h, 4,4) 



/ = 



-0.2500 0.2500 0.2500 

0.2500 0.7500 0.7500 

0.2500 0.7500 0.7500 

-0.2500 0.2500 0.2500 
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We may note that the Kernel of g 
i.e., {f\Af = 0} is given by 

V{I-S+S)z= f ZkVk 

k=K+\ 



Af, 



where z can be any arbitrary image. This can be used to 
obtain all the possible solutions of = Af by adding these 
arbitrary images to the GI solution. 

5. REGULARISATION 

We may easily note that A + A J and AA! + A J are no 
more singular if A > 0. Thus, we may compute: 

lambda= . 01 ; 

fh=inv (AtA+lambda^eye (size (AtA) ) ) ^ (A' ^p) ; 
reshape ( f h, 4,4) 



-0.2491 0.2497 0.2497 -0.2491 

0.2497 0.7484 0.7484 0.2497 

0.2497 0.7484 0.7484 0.2497 

-0.2491 0.2497 0.2497 -0.2491 



or still 

lambda= . 01 ; 

fh=A' ★inv (AAt + lambda^eye (size (AAt ) ) ) ★p; 
reshape ( f h, 4,4) 

which give the same solution. 

6. POSITIVITY CONSTRAINT 

We may remark that the problem is so ill-posed that, im- 
posing to the solutions to be of the minimum norme does 
not reduce enough the space of the possible solutions. The 
positivity constraint has been frequently used in many image 
reconstruction applications. A very simple technique to im- 
pose the positivity constraint in iterative algorithms is just to 
impose it at each iteration: 

for k=l:100 

fh=fh+. l^A' ^ (p-A^fh ( : ) ) ; 

fh=fh. ^ (fh>0) ; 

end 

reshape (fh, 4,4) ; 



fh. 



0.0000 0.0000 

0.0000 1.0000 1.0000 0.0000 

0.0000 1.0000 1.0000 0.0000 

0.0000 0.0000 



Of course, this technique is only one of the possible 
methods to use the prior information of the positivity. How- 
ever, we see that it can be very useful at least in this low 
scale case. But, as we will see later, in a real larger image 
reconstruction problem, it is not enough. 

7. REAL SIZE IMAGES IMPLEMENTATION OF 
THE ALGORITHMES 

Let consider a (256 x 256) pixel image. Then, it is no more 
question of really constructing the matrix A, because its di- 
mensions are (256^ x 256). Indeed, we do not really need its 
construction, we only need the results of the forward compu- 
tation Af and the backprojection A^g. The following pro- 
grams shows how to compute these quantities without actu- 
ally constructing the matrix A. 



function p=direct ( f ) ; 
pl=sum ( f ) ; 
p2=sum ( f ' ) ; 
p=[pl(:);p2(:)]; 
return 



function f=transp(p); 

l=length(p) ;pl=p (1:1/2) ; p2=p (1/2 + 1 : 1) ; 
f=ones (1/2, 1) ^pl' +p2^ones (1,1/2) ; 
return 

These programs can easily be used to obtain the follow- 
ing results/ 




a) UU b) 

Fig. 4 : a) Objet et ses projections, b) Resultat de 
retroprojection. 



We can also use them to impose any constraints such as 
positivity in the iterative methods: 

Moindre carre avec contraint de positivite : 



alpha= . 1 ; 

for k=l:100 

g=trans (p-direct (fh) ; 

f h=f h+alpha^g; 

fh=fh. ^ (fh>0) ; 

end 

Regularisation quadratique avec contraint de positivite ; 



alpha=.l;d=[-l -1;0 4 0;-l -1]; 
for k=l:100 

gO=trans (p-direct (fh) ; 
g=g0-lambda^conv2 (fh, d, ' same' ) ; 
f h=f h+alpha^g 
fh=fh. ^ (fh>0) ; 
end 

Remarquons que les algorithmes presentes plus haut 
sont assez rudimentaires (gradient a pas constant et a nom- 
bre d' iterations fini). II est evident que Ton peux faire 
mieux. A titre d'exemple, nous avons developpe un logi- 
ciel d' optimisation (gpave) un peu plus elabore qui met en 
oeuvre d' autre algorithmes d' optimisation, comme par exem- 
ple, gradient a pas adaptative, gradient conjugue et d'autres 
encore. 





a) 

Fig. 5 : a) Moindre carres avec contrainte de positivitee et 
b) Regularisation quadratique avec contrainte de positivitee. 



Les lignes de codes Matlab qui suivent montrent 1' usage 
de ce logiciel. II faut tout d'abord ecrire deux routines qui 
calculent le critere crit qui calcule 



reshape ( f h, 4,4) 



-M\Df\\ 



\\9-Af\ 

et son gradient dcrit 

V/ = -2A'{g - Af) + iXD'Df 

o\x Df correspond a 1' application d'une operation de con- 
volution de I'image /(/, j) avec une reponse impulsionnelle 

^ ^ ce qui correspond a 
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Avec ce programme d' optimisation il est alors facile de 
modifier les routines crit et dcrit pour changer le critere 
de la regularisation. Les figures suivantes montrent un cer- 
tain nombre des resultats. 



Notez aussi que D^Df correspond a 1' application d'une 
operation de convolution de I'image /(/, 7) avec une reponse 
impulsionnelle 



-1 -1 
0-4 
-1 -1 



ou * signifie une convolution. 

function J=crit ( fh, p, lambda) 

dp=p-direct (fh) ; 

JO = sum (dp ( : ) . 2 ) ; 

d=[-l 1;1 -1]; 

df=conv2 (fh, d, ' same' ) ; 

Jl=sum(df ( : ) .2) ; 

J= JO + lambda ^Jl ; 

return 

function dJ=dcrit ( fh, p, lambda) 

dp=p-direct (fh) ; 

dJ0=-2 ^transp (dp) ; 

d=[-l -1; 4 0;-l -1] ; 

dJl=conv2 (fh, d, ' same' ) ; 

dJ=dJO+lambda^dJl ; 

return 

Avec ces deux routines, le programme de la reconstruc- 
tion devient tres simple: 

f 0=transp (p) ; 
options = goptions; 
lambda=l ; 

fh=gpav ( ' crit' , f 0, options, 'dcrit' , p, lambda) 
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a b 
Fig. 6 : a) Moindre carres avec contraintes de positivitee et 
de support b) Regularisation quadratique avec contrainte de 
positivitee et de support. 



On remarque que le probleme est tres mal-conditionnee 
au sens que la manque d' information dans les donnees est 
trop important. II faut pouvoir obtenir d'autres donnees, 
i.e.,, des projections suivant d' autre angles. Pour cela 
il faut reecrire les routines directe et transp et les 
routines correspondantes crit et dcrit afin de pouvoir 
implementer le cas general du calcul des projections suivant 
n'importe quel angle et la retroprojection associee. 



Les figures suivantes montrent des exemples de recon- 
structions pour le cas ou on a 7 projections. 
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Fig. 7 : Reconstruction a partir de 7 projections: a) I'objet 
original et les donnees, b) retroprojection c) retroprojection 
filtre avec contrainte d) MC, e) MC avec positivite, f) MC 
avec positivite et contraint de support, g) Regularisation 
quadratique (RQ), h) RQ avec positivite, 



8. MODELISATION PAR CHAMPS DE MARKOV 
COMPOSITES 

Evidament, plus on a des donnees bien reparties, mieux 
sera les resultats. Mais, lorsque I'obtention d'autres pro- 
jections est impossible, il faudra recompenser la manque 
d' information par des modelisations plus precises. En parti- 
culier, dans le domaine du control non destructif (CND), une 
information a priori importante est que I'objet est compose 
d'un nombre fini de materiaux. Ceci signifie que 1' image que 
nous cherchons a reconstruire est composee d'un nombre fini 
de regions homogemes. C'est exactement la modelisation de 
cette information a priori qui est I'originalite des travaux que 
nous menons dans notre laboratoire. 

L'outil est la modelisation probabiliste par champs de 



Markov et 1' estimation bayesienne. Un grand nombre de 
travaux ont ete fait sur ce sujet (voir par exemple ifTOl fTTl 
fT2ll ). Ici, nous mentionnons seulement deux modelisations : 
Modelisation de 1' image par un champs composite (inten- 
sites-contours) ou (intensites-regions). Dans la premiere, 
on introduit une variable cachee binaire q{r) qui represent 
les contours et dans la deuxieme on introduit une variable 
cachee discrete z{r) qui peut prendre des valeurs discretes 
/: = 1 , • • • , /sT, representant les labels attribues aux pixels f{r) 
de r image ay ant les memes proprietes (par exemple se trou- 
vant dans une meme region homogene). 

8.1 Modele Intensites- Contours 

Dans cette modelisation, I'idee de base est de modeliser 
le fait qu'une image est en faite une fonction f{r) qui 
est continue par morceaux (piecewise continuous) ou par 
regions. II y a done des discontinuites (contours). On peut 
alors modeliser ces contours par une image binaire q{r). 
Le point essentiel est alors de decrire a I'aide d'une loi de 
probabilite conditionnelle p{f\q), le lien qu'il y entre des 
variables intensites / et des variables contours q qui peut 
etre resume par : 



Gas ID: 



Gas 2D: 



p{f{r)\q{r)J{s))=^\p{\-q{r)) 



3Gr(r) 



Ensuite, en choisissant une loi a priori appropriee pour p{q) 
et en utilisant des lois p{g\f) et p{f\q), on obtien la loi a 
posteriori p{f^q\g) qui peut etre utilisee pour inferer con- 
jointement sur / et sur q. A titre d' indication, considerons 
r estimation au sense du MAP : 

(/,^) =argmax {;?(/, gr|5f)} 

qui peut etre obtenu par un algorithme iteratif du type : 

/ = argmax/ {p{f\g, q)} = argminj {/(/)} 
q^=argmaxq {p{q\g)} 

avec 



J{f) = \\g-Hff^l^{l-q{r))[f{r)-l3 ^ f{s)] 

L'etape difficile est I'obtention de 1' expression de p{q\g) et 
surtout son optimisation, qui idealement ne peux se faire qu'a 
I'aide d'une recherche combinatoire. II existe un tres grand 
nombres de travaux portant sur differentes approximations 
qui permettent d'effectuer cette optimisation d'une maniere 
approchee mais realiste en cout de calcul pour des applica- 
tions reelles. 

Pour plus de detail sur cette methode se referer a ifTOlfTTl 
El. 

Ici, nous montrons un resultat typique que Ton peut 
obtenir avec de telles methode. Gomme nous pouvons con- 
state, cette modelisation a priori n'est pas encore suffisa- 
ment forte pour obtenir un resultat satisfaisant pour ce In- 
verse problem tres difficile. 
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a) UU b) 
Fig. 8 : Resultats de 1' estimation des intensites / en (a) et 
des contours q en (b). 



• MAP-Gibbs (Algorithm 2): 



8.2 Modele Intensites-Regions 

La modelisation precedente, bien que deja plus specifique, 
n'apportait pas d' information sur des valeurs qui peuvent 
prendre des pixels de 1' image. Dans certaines applications, 
par exemple en control non destructif (CND), nous savons a 
priori que I'objet est compose d'un nombre fini de materiaux 
(air, metal, composite). La modelisation qui suit permet de 
prendre en compte ce type d' information a priori. 

Plus specifiquement, nous proposons de modeliser 
r image par un champs composite (intensites-labels), ou les 
labels z{r) representent la nature de materiaux a la position 
du pixel r et I'homogeneite des intensites fk = {f{r)jr G 
^k} dans une region donnees = - z{r) = k} est 
modelise par un champs de Gauss-Markov: 

p{fk)=^{mkh^k) 
ce qui peut etre interprets aussi par des relations suivantes: 

p{f{r)\z{r)=k)=^{mk,Vk) 
P{f{r)) = Ef=i P{z{r) = k) ^{m,, a^) 
Eyt = diag[vi,--- ,vk\ 

ce qui montre que la distribution marginale de chaque pixel 
de r images est modelisee par un melange de gaussiennes. 

Supposant ensuite qu'a priori les pixels qui se trouvent 
dans deux regions differentes soient independantes, on peut 
ecrire : 

p{f\z) = f{^{mkl,l.k) 

k=\ 

La particularite de la methode que nous proposons est 
une modelization specifique pour des labels des regions 
(Champs de Potts) qui permet de decrire pour p{z) 



p{z) oc exp 



«I I mr)-z{s)) 



Notant aussi par = {cTe^, (m^^, O"^),/: = ^K} que 
Ton appelle le vecteur des hyperparametres, on aura a ex- 
primer la loi a posteriori jointe p{f,z^6\g), qui peut ensuite 
etre utilisee pour estimer ces inconnues. Differents choix 
sont alors possibles. Ici, nous mentionnons ceux que nous 
avons implementes et utilises: 
• MAP (Algorithm 1): 

/ =argmaxf{p{f\z,e,g)} 
e =^rgm^xe{p{0\f,z,g)} 
z = argmax;^ {p{z\ f, e,g)} 




= ^rgm^Xf{p{f\z,e,g)} 
avec p{e\f,z,g) 
avec p{z\f,0,g) 



Pour plus de details sur les expressions de ces lois et la mise 
en oeuvre de la methode dans un cadre plus general se referer 

adimisiiii. 

Principal avantage d'une telle modelisation et d'un tel 
methode est que Ton obtient non seulement une estimation 
de / mais aussi de z qui represente une segmentation de 
r image, et aussi par un simple algorithme de detection de 
contours sur z on obtiendrai aussi une image des contours q. 
La figure qui suit montre un resultat typique. 



c d 

Fig. 8 : a) Objet original et les deux projections, b,c et d) 
Resultats de 1' inversion par la methode proposee qui fourni 
une estimation des intensites (b), une estimation de labels en 
(c) et une estimation des contours en d). 



9. CONCLUSION 

Dans ce travail, a but pedagogique, au travers d'un Inverse 
problem de la reconstruction d' image en Tomography X 
lorsque le nombre de projections sont tres limite, nous avons 
analyse les difficultes inherentes des problemes inverses. 
Le principal objectif etait de montrer que les differentes 
methodes classiques naives, mais tres utilisees, ne donnent 
pas de solutions satisfaisantes et qu'il y a un besoin de pro- 
poser des methodes d' inversion plus sophistiquees qui per- 
mettent d'introduire de 1' information a priori necessaire pour 
compenser la manque d' information dans les donnees. 

Un grand nombre de modelisations ont ete proposees 
(voir par exemple ifTTl fTSl ). Mais, ici, nous nous sommes 
contente des methodes qui modelisent 1' image au niveau des 
pixels. 

Le cas particulier de la reconstruction a partir de 
deux projections est detaille et une methode basee sur la 
modelisation de 1' image par un champs de Markov compose 
(intensite-labels) et 1' estimation bayesienne est presentee qui 



permet, au moins, d'obtenir une solution satisfaisante au 
pro bleme. Les routines Matlab correspondant est disponible 
sur |http : //djaf ari . fre e . f r/TomoX[ 
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ABSTRACT 

Ce travail, a but pedagogique, presente le probleme inverse 
de la reconstruction d' image en tomographic X lorsque le 
nombre des projections est tres limite. L'objectif est de mon- 
trer que les differentes methodes classiques naives, mais tres 
utilisees pour sa resolution, ne donnent pas de solutions sat- 
isfaisantes. II est alors necessaire de proposer des methodes 
d' inversion plus sophistiquees qui permettent d'introduire de 
r information a priori necessaire pour I'obtention d'une so- 
lution acceptable. Le cas particulier de la reconstruction 
a partir de deux projections est detaille, les resultats que 
Ton peux obtenir avec les differentes methodes algebriques 
sont presentes et la necessite des methodes probabilistes 
est demontree. Finalement, une methode basee sur la 
modelisation de 1' image par un champs de Markov compose 
(intensite-labels) et 1' estimation bayesienne est presentee qui 
permet, au moins, d' obtenir une solution satisfaisante au 
probleme. Dans ce travail, 1' accent est mise sur 1' aspect 
pedagogique. En effet, ce probleme est traite dans le cadre 
d'un enseignement au niveau du DEA et les outils (pro- 
grammes Matlab) necessaires pour une demonstration seront 
disponibles lors de 1' expose. 



1. INTRODUCTION 

Le probleme de la reconstruction d' image en tomographic X 
est presente dans Fig. 1 et Fig. 2 montre la modelisation du 
probleme via la Transformee de Radon et le schema de la 
methode classique de retroprojection filtree d [2l [3l (H [S] [6l 
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(n,r2) = / f{x,y,z) dl g^{r) = / f{x,y) dl 



Probleme directe: f{x^y) or f{x,y^z) 
Probleme inverse: (r) or (ri , r2) 



' f{x,y) oxf{x,y,z) 



Fig. 1 : Tomographie X 




Derivation 



g(r, 0) = jj f{x^y)8{r — xco^(^—y^m(^)dxdy 
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Fig. 2 : Tomographie X, Transformee de Radon et 
Methode analytique de Retroprojection filtree 



II est evident que ces methodes analytiques 
(Retroprojection ou Retroprojection filtree) ne donne 
des resultats satisfaisants que lorsqu'il y a un grand nom- 
bre de projections. C'est pourquoi, dans ce travail, nous 
considerons les methodes algebriques qui permettent plus 
de souplesse pour le developpement des methodes plus 
sophistiquees. Fig. 3 montre I'etape de la discretisation du 
probleme qui le transforme a celui de la resolution d'un 
systeme d' equations lineaires g = Hf. Nous montrons 
ensuite que le probleme est sous-determine et qu'il y a 
une infinite de solutions possibles. Plus interessant encore 
qu'aucune methode algebrique du type Moindres carres 
(MC), Inversion generalisee, ou meme la regularisation 
quadratique E [51, ne fournit une solution satisfaisante. 
Aussi, r application des contraintes de positivite, bien 



qu'ameliore les resultats de ces methodes, ne suffit pas, d'ou 
la necessite de proposer des modelisations a priori plus 
complexes. 

2. DISCRETIZATION DU PROBLEME 

Dans ce travail, afin de pedagogie, nous considerons un 
probleme de reconstruction d' image avec de dimensions 
reduites a partir de deux projections horizontale et verti- 
cale, ce qui permet de mieux apprehender les difficultes du 
probleme. 
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Fig. 3 : Discretisation du probleme 



Soient / = [/i , • • • ,/i6]^ les valeurs des pixels de I'image 
f{x^y) de dimensions (4 x 4) et = [^i, • • • ,^3]^ l^s valeurs 
de ses projections g{r^ (^>) suivant les angles = et = 90 
degres. Supposons Ax = 1, A_y = 1, Ar = 1 et considerons les 
deux representations suivantes: 
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Notons aussi 
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g4f = fell,- 
>gsf = [g2U- 

et formons les matrices Ai , A2 et A telle qu'on puisse ecrire 



gi=Aif, g2 
Considerons I'image 



A2f, g = Af 
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et calculons sa projection en utilisant les lignes Matlab suiv- 
antes : 

Al=[ 1111000000000000; 
0000111100000000; 
0000000011110000; 
0000000000001111]; 

A2=[ 0001000100010001; 
0010001000100010; 
0100010001000100; 
1000100010001000]; 

A= [Al; A2] ; 

f = [ ; 1 1 ; 1 1 ; ] ; 

P=A^f (:) ; 
ce qui donne : 

g' = [02200220] 

Ainsi, la resolution du probleme direct ne pose aucune dif- 
ficulte. ]V[ais, considerons maintenant le probleme inverse: 
etant donnee g trouver /. 

3. INDETERMINATION DU PROBLEIME 

II est evident que ce probleme a une infinite de solutions pos- 
sibles. En voici quatre: 
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-.5 .5 
12 0-1 
-10 2 1 
0.5 -.5 





2 

2 





-.5 .5 

2 

2 
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4. EQUIVALENCE ALGEBRIQUE DE 
RETROPROJECTION 

En comparant les relation continue et discrete et les 
operateurs direct et adjoint, on trouve que 1' operation de 
retroprojection en continue correspond a I'operateur de trans- 
position de matrice en discret. Ainsi, la solution / = A^g 
correspond a la solution au sense de la la retro-projection. 



4^ 
^2 



et calculons 



Exprimons alors la matrice = [ A\ 
cette solution : 

fh=A' ★p; reshape (fh, 4, 4) 

2 2 

2 4 4 2 

2 4 4 2 

2 2 



Notons que / = A^g = A[gi + ^25^2 est 1' addition de deux 
images 



2 2 

2 2 

2 2 

2 2 





2 2 2 2 

2 2 2 2 





chacune etant la retro-projection d'une des deux projections. 



Remarquons aussi que cette image, a une constant pres 
est tres proche du resultat de la convolution de 1' image 
d'origine avec la reponse impulsionnelle 



1 

1 2 1 
1 



6. MOINDRES CARRES DE NORME MINIMAL, 
INVERSION GENERALISEE ET DTVS 

Notons que ces deux matrices sont singulieres. Rappelons 
qu'une solution au sens des moindres carres s'ecrit: 

/ = argmin{||g-A/||2}, 



5. INVERSION GENERALISEE 

En vue de la definition d'une solution au sens d' inversion 
generalisee, calculous les matrices A' A et AA': 



AA' = 
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4000 1111 
0400 1111 
0040 1 1 1 1 
0004 1 1 1 1 
1 1 1 1 4000 
11110400 
1 1 1 1 0040 
1 1 1 1 0004 
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A',A2 _ 



A[Ai 



1 + / 



A[A2 
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2 111 
12 11 
112 1 
1112 
1000 
0100 
00 1 
000 1 
1000 
0100 
00 1 
000 1 
1000 
0100 
0010 
000 1 



1000 
100 
00 1 
000 1 
2 111 
12 11 
112 1 
1112 
1000 
100 
00 1 
000 1 
1000 
100 
00 1 
000 1 



1000 
100 
00 1 
000 1 
1000 
100 
00 1 
000 1 
2 111 
12 11 
112 1 
1112 
1000 
100 
00 1 
000 1 



1000 
100 
00 1 
000 1 
1000 
100 
00 1 
000 1 
1000 
100 
00 1 
000 1 
2 111 
12 11 
112 1 
1112 



Calculons les valeurs singulieres des matrices AA^ A^ A: 

AAt=A^A' ; svd ( AAt ) 



et si la matrice A^A est inversible on obtient : 
{A'A)-^A'g. 

De meme, une solution de norme minimale est 

/ = arg min {\\ff} 

Af=g 



et si la matrice AA^ est inversible on obtient : / = 
A'{AA')-^g. 

Nous ne pouvons alors calculer aucune de ces deux 
solutions car aucune des deux matrices AA^ et A^A est 
inversible. Notons cependant que si on ne garde que les 
elements diagonaux de ces deux matrices, on obtient des 
resultats suivants : 

fh=diag (1 . /diag (At A) ) ★A' ^p; reshape (fh, 4, 4) 

110 
12 2 1 
12 2 1 
110 

fh=A' ^diag (1 . /diag (AAt) ) ★p; reshape (fh, 4,4) 

.5 .5 
.5 1 1 .5 
.5 1 1 .5 
.5 .5 

II est cependant possible de calculer la solution inverse 
generalisee qui est la solution de norme minimale de 
A^Af = A^g en utilisant la decomposition tronquee des 
valeurs singulieres (DTVS): 



/ = 



/ = L — ^ — 



k=l 



ou Uk et Vk sont, respectivement, des vecteurs propres de 
AA^ et de A^A et sont des valeurs singulieres associees. 

[U, S,V]=svd(A) ; 

s=diag (S) ; s 1= [ 1 . / s ( 1 : 7 ) ; zeros (1,1)]; 
Sl= [diag (si) ; zeros (8,8)]; 
fh=V^Sl^U' ^p; reshape (fh, 4, 4) 

Dans cet exemple ^ = 7 et la solution IG peut etre cal- 
culee par : 

f h=svdpca (A, p, . 1 , 7 ) ; reshape ( f h, 4,4) 



svd(AAO = [84444440] 

AtA=A' ^A; svd (AtA) ; 

svd(A'A) = [8444444000000000] 



-0.2500 0.2500 0.2500 -0.2500 

0.2500 0.7500 0.7500 0.2500 

0.2500 0.7500 0.7500 0.2500 

-0.2500 0.2500 0.2500 -0.2500 



ou encore par ralgorithme iteratif suivant: 

for k=l:100; 

fh=fh+. l^A' ^ (p-A^fh ( : ) ) ; 
end; 

reshape ( f h, 4,4) 



-0.2500 0.2500 0.2500 -0.2500 

0.2500 0.7500 0.7500 0.2500 

0.2500 0.7500 0.7500 0.2500 

-0.2500 0.2500 0.2500 -0.2500 



Notons aussi que le noyau de la transformation lineaire 
g = AfJ.e., {/|A/ = 0}est 



V{I-S^S)z= ^ ZkVk 

k=K+l 

avec z un vecteur arbitraire. Ceci nous permet de trouver 
toutes les solutions possibles du probleme en rajoutant ce 
terme arbitraire a la solutions IG. 

7. REGULARISATION 

Notons que, A + A/ et AA^ + A J sont inversibles pour 
A > 0. Ceci nous permet de calculer 

lambda= . 01 ; 

fh=inv (AtA+lambda^eye (size (At A) ) ) ^ (A' ^p) ; 
reshape ( f h, 4,4) 



-0.2491 0.2497 0.2497 -0.2491 

0.2497 0.7484 0.7484 0.2497 

0.2497 0.7484 0.7484 0.2497 

-0.2491 0.2497 0.2497 -0.2491 



ou encore 



lambda= . 01 ; 

fh=A' ^inv (AAt + lambda^eye (size (AAt ) ) ) ^p; 
reshape ( f h, 4,4) 

qui fourni la meme solution. 

8. CONTRAINTE DE POSITIVITE 

On peut remarquer que la manque d' information dans les 
donnees est telle que la contrainte de norme minimal ne 
restreint pas suffisament I'espace des solutions possibles. 
Dans les problemes inverses en imagerie, une information 
qui est souvent disponible est la positivite de la solution. 
Imposer alors a la solution d'etre positive est alors une 
technique souvent utilisee. Une approche simple dans les 
methodes iterative pour imposer cette contraintes est simple- 
ment 1' imposer a chaque iteration : 

for k=l:100 

fh=fh+. l^A' ^ (p-A^fh ( : ) ) ; 

fh=fh. ^ (fh>0) ; 

end 

reshape (fh, 4,4) ; 



fh 





0.0000 
0.0000 




0.0000 0.0000 
1.0000 1.0000 0.0000 
1.0000 1.0000 0.0000 



0.0000 0.0000 







Bien entendu, ceci n'est qu'une methode simple et il ex- 
iste un grand nombre d'algorithmes d' optimisation sous con- 
traintes que Ton peut utiliser, mais la description de ces al- 
gorithmes sort du cadre de ce travail. 

9. MISE EN OEUVRE DANS UN CAS REEL 

Examinons maintenant le cas d'une image de plus grande 
taille (256 x 256). Ici, il n'est pas question de former la 
matrice A car de dimensions (256^ x 256). Par contre, on 
aura besoin de calculer Af et A^g, mais pour cela on n'a 
pas besoins de construire reellement la matrice A. Les deux 
fonctions suivantes effectuent ces deux taches : 

function p=direct(f); 
pl=sum ( f ) ; 
p2=sum ( f ' ) ; 
p=[pl(:);p2(:)]; 
return 

function f=transp(p); 

l=length(p) ;pl=p (1:1/2) ; p2=p (1/2 + 1 : 1) ; 
f=ones (1/2,1) ^pl' +p2^ones (1,1/2); 
return 

Ces routines peuvent alors facilement utilisees pour 
obtenir les images de la figure qui suit. 



a) 




Fig. 4 : a) Objet et ses projections, b) Resultat de 
retroprojection. 



Nous pouvons aussi utiliser ces routines pour mettre en 
oeuvre la pluparts des methodes iteratives : 

Moindre carre avec contraint de positivite : 

alpha= . 1 ; 
for k=l:100 
g=trans (p-direct (fh) ; 
f h=f h+alpha★g; 
fh=fh. ^ (fh>0) ; 
end 

Regularisation quadratique avec contraint de positivite : 

alpha=.l;d=[-l -1; 4 0;-l -1]; 
for k=l:100 



gO=trans (p-direct (fh) ; 
g=g0-lambda^conv2 (fh, d, ' same' ) ; 
f h=f h+alpha^g 
fh=fh. ^ (fh>0) ; 
end 

Remarquons que les algorithmes presentes plus haut 
sont assez rudimentaires (gradient a pas constant et a nom- 
bre d' iterations fini). II est evident que Ton peux faire 
mieux. A titre d'exemple, nous avons developpe un logi- 
ciel d' optimisation (gpave) un peu plus elabore qui met en 
oeuvre d' autre algorithmes d' optimisation, comme par exem- 
ple, gradient a pas adaptative, gradient conjugue et d'autres 
encore. 




Fig. 5 : a) Moindre carres avec contrainte de positivitee 
et b) Regularisation quadratique avec contrainte de 
positivitee. 



Les lignes de codes Matlab qui suivent montrent 1' usage 
de ce logiciel. II faut tout d'abord ecrire deux routines qui 
calculent le critere crit qui calcule 

J=\\g-Aff + ^Dff 

et son gradient dcrit 

V/ = -2A'{g - Af) + iXD'Df 

ou Df correspond a 1' application d'une operation de convo- 
lution de I'image f{ij) avec une reponse impulsionnelle 

" -1 1 " 
1 -1 

ce qui correspond a 

EE -/('■- hj)? + \f{.iJ)-f{.iJ-\)?) ■ 

i j 

Notez aussi que D^Df correspond a 1' application d'une 
operation de convolution de I'image /(/, 7) avec une reponse 
impulsionnelle 



ou * signifie une convolution. 

function J=crit ( fh, p, lambda) 
dp=p-direct (fh) ; 
JO=sum(dp ( : ) .2) ; 



d=[-l 1;1 -1]; 
df=conv2 (fh, d, ' same' ) ; 
Jl=sum(df ( : ) .2) ; 
J= JO+lambda^ Jl ; 
return 

function dJ=dcrit ( fh, p, lambda) 

dp=p-direct (fh) ; 

dJ0=-2 ★transp (dp) ; 

d=[-l -1;0 4 0;-l -1] ; 

dJl=conv2 (fh, d, ' same' ) ; 

dJ=dJO+lambda^dJl ; 

return 

Avec ces deux routines, le programme de la reconstruction 
devient tres simple: 

f 0=transp (p) ; 
options = goptions; 
lambda=l ; 

fh=gpav ( ' crit' , f 0, options, 'dcrit' , p, lambda) ; 
reshape ( f h, 4,4) 

Avec ce programme d' optimisation il est alors facile de 
modifier les routines crit et dcrit pour changer le critere 
de la regularisation. Les figures suivantes montrent un cer- 
tain nombre des resultats. 




a b 



Fig. 6 : a) Moindre carres avec contraintes de positivitee 
et de support b) Regularisation quadratique avec 
contrainte de positivitee et de support. 



On remarque que le probleme est tres mal-conditionnee 
au sens que la manque d' information dans les donnees est 
trop important. II faut pouvoir obtenir d'autres donnees, 
i.e.,, des projections suivant d' autre angles. Pour cela 
il faut reecrire les routines directe et transp et les 
routines correspondantes crit et dcrit afin de pouvoir 
implementer le cas general du calcul des projections suiv- 
ant n'importe quel angle et la retroprojection associee. Ces 
routines sont bien sure plus complexes car elles ont besoin 
de toutes les parametres geometriques; les positions et les 
nombre des sources et des capteurs, les positions et les pas 
de discretisation suivant les axes de I'objets, etc. 

Nous avons developpe un ensemble de routines Matlab, 
sous la forme d'un toolbox [?], qui permet de simuler les situ- 
ations de geometric paralelle ou conique pour des operateurs 
de projections et retroprojection, ainsi que la mise en oeuvre 
des methodes de reconstructions deja mentionnees dans ce 
contexte general. Les figures suivantes montrent des exem- 
ples de reconstructions pour le cas ou on a 7 projections. 




Fig. 7 : Reconstruction a partir de 7 projections: a) 
I'objet original et les donnees, b) retroprojection c) 
retroprojection filtre avec contrainte d) MC, e) MC avec 
positivite, f) MC avec positivite et contraint de support, 
g) Regularisation quadratique (RQ), h) RQ avec 
positivite. 



10. MODELISATION PAR CHAMPS DE MARKOV 
COMPOSITES 

Evidament, plus on a des donnees bien reparties, mieux 
sera les resultats. Mais, lorsque I'obtention d'autres pro- 
jections est impossible, il faudra recompenser la manque 
d' information par des modelisations plus precises. En parti- 
culier, dans le domaine du control non destructif (CND), une 
information a priori importante est que I'objet est compose 
d'un nombre fini de materiaux. Ceci signifie que 1' image que 
nous cherchons a reconstruire est composee d'un nombre fini 
de regions homogemes. C'est exactement la modelisation de 
cette information a priori qui est I'originalite des travaux que 
nous menons dans notre laboratoire. 



L'outil est la modelisation probabiliste par champs de 
Markov et 1' estimation bayesienne. Un grand nombre de 
travaux ont ete fait sur ce sujet (voir par exemple ifTOl fTTl 
fT2l ). Ici, nous mentionnons seulement deux modelisations : 
Modelisation de 1' image par un champs composite (inten- 
sites-contours) ou (intensites-regions). Dans la premiere, 
on introduit une variable cachee binaire q{r) qui represent 
les contours et dans la deuxieme on introduit une variable 
cachee discrete z{r) qui peut prendre des valeurs discretes 
/: = 1 , • • • , /sT, representant les labels attribues aux pixels /(r ) 
de r image ay ant les memes proprietes (par exemple se trou- 
vant dans une meme region homogene). 

10.1 Modele Intensites- Contours 

Dans cette modelisation, I'idee de base est de modeliser 
le fait qu'une image est en faite une fonction f{r) qui 
est continue par morceaux (piecewise continuous) ou par 
regions. II y a done des discontinuites (contours). On peut 
alors modeliser ces contours par une image binaire q{r). 
Le point essentiel est alors de decrire a I'aide d'une loi de 
probabilite conditionnelle p{f\q), le lien qu'il y entre des 
variables intensites / et des variables contours q qui peut 
etre resume par : 

Gas ID: 

p{fj\qjJiJ7^j)=^{fj\P{l-qj)fj-U0j) 

OU nous avons utilise la notation ^(x|m, v) pour representer 
une distribution gaussienne de la variable x avec la moyenne 
m et la variance v. 

Gas 2D: 

p{f{rMr)J{s))=^(f{r)\P{l-q{r)) ^ f{s),aj 

Ensuite, en choisissant une loi a priori appropriee pour p{q) 
et en utilisant des lois p{g\f) et p{f\q), on obtien la loi a 
posteriori p{f,q\g) qui peut etre utilisee pour inferer con- 
jointement sur / et sur q. A titre d' indication, considerons 
r estimation au sense du MAP : 

(/,^) =argmax {/?(/, gr|gf)} 

qui peut etre obtenu par un algorithme iteratif du type : 

/ = argmaxj {p{f\g, q)} = argminj {/(/)} 
q^=argmaxq {p{q\g)} 

avec 

J{f) = \\g-Hff + l^{l-q{r))(f{r)-p ^ f{s)] 

L'etape difficile est I'obtention de 1' expression de p{q\g) et 
surtout son optimisation, qui idealement ne peux se faire qu'a 
I'aide d'une recherche combinatoire. II existe un tres grand 
nombres de travaux portant sur differentes approximations 
qui permettent d'effectuer cette optimisation d'une maniere 
approchee mais realiste en cout de calcul pour des applica- 
tions reelles. 



Pour plus de detail sur cette methode se ref erer a |[TOl [TTJ 
H. 

Ici, nous montrons un resultat typique que Ton peut 
obtenir avec de telles methode. Comme nous pouvons con- 
stater, cette modelisation a priori n'est pas encore suff- 
isament forte pour obtenir un resultat satisfaisant pour ce 
probleme inverse tres difficile. 
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Fig. 8 : Resultats de restimation des intensites / en (a) et 
des contours q en (b). 



10.2 Modele Intensites-Regions 

La modelisation precedente, bien que deja plus specifique, 
n'apportait pas d' information sur des valeurs qui peuvent 
prendre des pixels de 1' image. Dans certaines applications, 
par exemple en control non destructif (CND), nous savons a 
priori que I'objet est compose d'un nombre fini de materiaux 
(air, metal, composite). La modelisation qui suit permet de 
prendre en compte ce type d' information a priori. 

Plus specifiquement, nous proposons de modeliser 
r image par un champs composite (intensites-labels), ou les 
labels z{r) representent la nature de materiaux a la position 
du pixel r et I'homogeneite des intensites fk = {f{r).,r G 
^k] dans une region donnees = • z{r) = k} est 
modelise par un champs de Gauss-Markov: 

p{fk) = ^{f\mkl,l.k) 
ce qui peut etre interprets aussi par des relations suivantes: 

p{f{r)\z{r)=k)=^{f{r)\mk,Vk) 

P{f{r))=n- 
= diag [vi,- 



p{f{r))=l^Up{z{r)=k) ^{f{r)\m,,Gl) 



ce qui montre que la distribution marginale des pixels de 
r images est modelisee par un melange de gaussiennes. 

Supposant ensuite qu'a priori les pixels qui se trouvent 
dans deux regions differentes soient independantes. On peut 
alors ecrire : 

p{f\z) =nf=i^(/Ki,^^) 

= \\l=i\\reR,^{f[r)\mk,Vk) 

La particularite de la methode que nous proposons est 
une modelization specifique pour des labels des regions 
(Champs de Potts) qui permet de decrire pour p{z) 



p{z) oc exp 



a 



-z{s)) 



Notons par Q = {cjg^, (m^^, ol)^k = 1, • • • ,^}, appelle le 
vecteur des hyperparametres. Dans un probleme reel (recon- 
struction non supervise), ce vecteur est aussi inconnu et il 
faut I'estimer. II faut alors lui attribuer une loi a priori p{0). 
Une fois fait, nous avons tous les elements en main pour ex- 
primer la loi a posteriori joints 

p{f,z,e\g) ocp{f\z,e,g)p{z)p{e) 

o^p{g\f.e)p{f\z,e)p{z)p{o). 

On peut ensuite utiliser cette loi pour inferer toutes ces in- 
connues. Cette etape se fait en general par I'intermediaire 
de la definition des estimateurs /, £ et 6 qui correspondent 
soit au mode ou la moyenne de cette lois a posteriori, qui 
ne peuvent, en general etre calcules que par des algorithmes 
iteratives qui utilisent successivement les lois a posteriori 
conditionnelles 



p{f\z,d,g)<^ 
p{z\f,d,g)<^ 
p{e\f,z,g)c< 



p{9\f,e)p{f\z,d) 

p{z\f,d)p{z) 

p{f\z,e,g)p{e) 



Differents choix sont alors possibles. Ici, nous mention- 
nons ceux que nous avons implementes et utilises: 

• MAP (Algorithm 1): 

: argmax/{/?(/|2;,e,5f)} 
: ^Ygmaxz{p{z\f,e,g)} 
-- ^rgm^xe{p{e\f,z,g)} 



• Gibbs (Algorithm 2): 

echant. / 

echant. z 

echant. 6 



avec p{f\z,e,g) 
avec p{z\f,e,g) 
avec p{0\f,z,g) 



MAP-Gibbs (Algorithm 3): 

/= argmaxj{;?(/|2;,e,5f)} 
echant. z avec p{z\f,6,g) 
echant. avec p{0\f^z^g) 

Gibbs-EM (Algorithm 4): 

/= E{f\z,e,g} = ^rgm^Xf{p{f\z,e,g)} 
z= E{z\f,e,g} ou echant. avec /?(2; I/, e,5f) 
0= E{0\f^z^g} ou echant. avec p{0\f^z^g) 

L' equivalence de la premiere ligne est due au fait que la lois 
a posteriori p{f\z^6^g)est une gaussienne et done sa mode 
et sa moyenne se confondent. 

Pour plus de details sur les expressions des lois condi- 
tionnelles qui interviennent dans ces algorithmes et la mise 
en oeuvre de la methode dans un cadre plus general se referer 

aiiniiiiEiiisi. 

Principal avantage d'une telle modelisation et d'un tel 
methode est que Ton obtient non seulement une estimation 
de / mais aussi de z qui represente une segmentation de 
r image, et aussi par un simple algorithme de detection de 
contours sur z on obtiendrai aussi une image des contours q. 
La figure qui suit montre un resultat typique. 
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Fig. 9 : a) Objet original et les deux projections, b,c et d) 
Resultats de I'inversion par la methode proposee qui 
fourni une estimation des intensites (b), une estimation 
de labels en (c) et une estimation des contours en d). 



11. CONCLUSION 

Dans ce travail, a but pedagogique, au travers d'un probleme 
inverse de la reconstruction d' image en tomographic X 
lorsque le nombre de projections sont tres limite, nous avons 
analyse les difficultes inherentes des problemes inverses. 
Le principal objectif etait de montrer que les differentes 
methodes classiques naives, mais tres utilisees, ne donnent 
pas de solutions satisfaisantes et qu'il y a un besoin de pro- 
poser des methodes d' inversion plus sophistiquees qui per- 
mettent d'introduire de 1' information a priori necessaire pour 
compenser la manque d' information dans les donnees. 

Un grand nombre de modelisations ont ete proposees 
(voir par exemple |[T7| [El). Mais, ici, nous nous sommes 
contente des methodes qui modelisent 1' image au niveau des 
pixels. 

Le cas particulier de la reconstruction a partir de 
deux projections est detaille et une methode basee sur la 
modelisation de 1' image par un champs de Markov compose 
(intensite-labels) et 1' estimation bayesienne est presentee qui 
permet, au moins, d'obtenir une solution satisfaisante au 
pro bleme. Les routines Matlab correspondant est disponible 
sur |http ://djafari.free. f r/TomoX| 
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